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ABSTRACT 


A technique  is  presented  that  allows  the  numerical  evaluation  of 
the  velocity  potential  for  a pulsating  source  in  a fluid  with  a free 
surface  and  finite  depth.  The  series  expansion  used  to  obtain  a 
numerically  tractable  representation  of  the  finite  depth  potential  is 
in  a form  suitable  for  use  with  ship  motion  programs  that  employ  source 
distribution  techniques  in  conjunction  with  strip  theory.  The  effects 
of  water  depth,  source  submergence,  and  pulsating  frequency  on  the 
free-surface  elevation  are  explored  with  the  intent  of  demonstrating 
finite  depth  effects  on  a readily  identifiable  physical  parameter. 

ADMINISTRATIVE  INFORMATION 

The  study  described  herein  evolved  from  research  supported  by  the  Naval 
Ship  Systems  Command  (NAVSHIPS)  under  the  General  Hydromechanics  Research 
Program  of  the  Naval  Ship  Research  and  Development  Center  (NSRDC) . Funding 
was  provided  under  Subproject  S-R02301,  Task  OlOl. 


INTRODUCTION 


One  of  the  most  useful  techniques  employed  In  the  theoretical  analysis 
of  ship  motions  Is  the  "strip  theory"  approach,  originally  derived  for  pitch 
and  heave  motions  by  Korvln-Kroukovsky  and  W.R.  Jacobs  (1957)  and  later 
improved  and  extended  by  Salvesen,  Tuck,  and  Faltinsen  (1970).  This  theory  is 
based  on  the  assumption  that  the  three-dimensional  hydrodynamic  properties  of 
a ship  may  be  expressed  as  the  sum  of  the  two-dimensional  characteristics  of 
specified  transverse  cross  sections  of  the  ship  hull.  Strip  theory  has 
proved  to  be  a useful  tool  for  analytic  ship  motion  analysis  and  enjoys  wide 
application.  Since  the  use  of  the  theory  requires  that  the  two-dimensional 
nydrodynamic  characteristics  of  various  cross-sections  be  determined, 
techniques  designed  to  provide  this  Information  form  an  important  role  for 
proper  application  of  strip  theory  to  ship  motion  work.  The  two  principal 
means  of  obtaining  these  properties  are  the  multipole  expansion  method 
(often  referred  to  as  the  conformal  mapping  technique)  and  the  source 
distribution  method. 

The  multipole  expansion  technique  was  first  employed  by  Ursell  (1949) 
for  determining  the  hydrodynamic  characteristics  of  a circular  cylinder 
oscillating  in  the  free  surface.  Several  investigators  have  contributed  to 
the  extension  of  this  technique  to  the  Lewis  forms  and  later  to  more  general 
shapes;  especially  notable  are  the  works  of  Grimm  (1953)  and  Porter  (1960). 
Korvin-Kroukovsky  (1957),  Gerrltama  (1966),  Smith  (1966),  Grim  (1960),  and 
Vassllopoulos  (1964)  all  employed  the  multipole  expansion  technique  to 
determine  the  two-dimensional  section  characteristics  for  use  with  strip 
theory . 
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The  source  distribution  technique  was  developed  by  Frank  (1967).  The 
Frank  approach  has  been  very  successful  because  of  the  wide  variety  of 
unusual  forms  that  can  be  treated,  e.q.  bulbous  bows,  submerged  bodies, 
and  catamarans.  Thus,  certain  bulbous  bow  shapes  that  are  difficult  to 
accurately  represent  with  mapping  techniques  offer  no  special  difficulties 
for  the  close-fit  source  distribution  approach.  The  source  distribution 
technique  has  been  used  in  the  most  advanced  ship  motion  work  currently 
available.  The  work  of  Frank  and  Salvesen  (1970)  utilizes  this  technique 
for  pitch  and  heave  motion  prediction.  The  most  recent  effort  in  the 
field  by  Salvesen,  Tuck,  and  Faltinsen  (1970)  employs  the  source  distribution 
technique  to  determine  the  hydrodynamic  characteristics  of  a ship  in  all 
six  degrees  of  freedom. 

The  original  work  by  Frank  is  restricted  to  the  deep-water  case. 

However,  problems  such  as  operation  of  supertankers  within  near-shore 
waterways  and  naval  vessels  close  to  shore  have  created  interest  in  the 
behavior  of  ships  in  finite  depth  waters.  It  is  desirable,  therefore, 
to  extend  the  infinite  depth  fluid  case  to  include  the  effects  of  water 
of  finite  depth.  The  first  and  most  difficult  requirement  for  this 
extension  is  to  obtain  a numerical  method  for  predicting  the  velocity 
potential  for  a pulsating  source  of  arbitrary  strength  in  a fluid  with 
a free  surface  and  finite  depth.  With  a suitable  evaluation  of  this 
potential,  the  next  step  is  to  adapt  it  for  use  with  a source  distribution 
program  such  as  that  by  Frank  (1967).  This  will  allow  two-dimensional 
added-mass  and  damping  characteristics  to  be  predicted  for  a variety  of 
shapes.  Such  data  can  then  be  used  with  strip  theory  to  determine 
motions  of  various  marine  vehicles  in  a finite  water  depth  environment. 
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The  objective  of  this  work  is  the  development  of  a numerical  technique 
to  accomplish  the  first  of  the  ebove  steps.  Specifically,  the  following 
efforts  are  presented: 

1.  Derivation  of  the  velocity  potential  for  a single  pulsating 

source  in  a fluid  with  a free  surface  and  finite  depth. 

2.  A numerical  technique  for  evaluation  of  the  resultant  potential. 

3.  Exploration  of  the  following  parameters: 

a.  Water  depth  h 

b.  Source  submergence  b 

c.  Pulsating  frequency  o 

In  addition,  it  will  be  demonstrated  that  as  the  water  depth  Increases,  the 
solution  tends  toward  that  obtained  for  infinitely  deep  fluids. 

Ths  required  finite  depth  velocity  potential  has  been  given  by  several 
investigators:  John  (1950),  Wehausen  and  Laltone  (1960),  Thorne  (1953)  and 
Porter  (1960).  Eoth  Wehausen  and  Laltone  (1960)  and  also  Porter  (1960) 
provide  detailed  explanations  of  the  procedure  used  to  formulate  the  velocity 
potential  for  both  the  Infinite  and  finite  depth  two-dimensional  problems. 
Thorne  (1953)  also  provides  a good  description  of  the  procedures  used  to 
formulate  the  potential  but  gives  considerably  less  detail. 

Numerical  evaluations  of  the  finite  depth  potential  for  the  purpose  of 
finding  added  mass  and  damping  have  been  attempted  by  several  investigators 
with  mixed  success.  Yu  and  Ursell  (1961)  used  infinite  series  expansions  of 
a type  suggested  by  Thome  (1953)  to  allow  comparison  of  experimental  and 
theoretical  wave  amplitudes  and  obtained  generally  good  agreement.  McCreight 
(1970)  experienced  some  difficulty  with  his  numerical  evaluation  as  depth 
Increased,  but  he  obtained  good  results  for  shallow  water  values  of  his  depth 
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parameters.  Portar  (1960)  axparlencad  difficulty  In  the  numerical  evalu- 
ation of  hie  finite  depth  potential  and  had  to  abandon  the  effort  before 
the  problame  could  be  remedied.  Kim  (1969)  appeared  to  ha/e  some  diffi- 
culty in  comparing  hia  result*  with  those  of  Yu  and  Ursell  (1961),  and  his 
effort  may  contain  soma  errors.  Kim  states  that  he  employed  the  same 
expansions  of  the  velocity  potential  as  used  by  Yu  and  Ursell  and  one 
would  expect  their  results  to  be  in  close  agreement. 

The  technique  used  in  this  report  to  evaluate  the  velocity  potential 
follows  an  approach  due  primarily  to  John  (1950)  and  \c.v_u3~n  i-'.d  ltene 
(1960),  although  the  derivation  of  the  potential  in  Appendix  A follows  an 
approach  suggested  by  Thorne  (1953).  The  form  of  the  potential  ultimately 
employed  for  numerical  evaluation  is  slightly  different  than  that  suggested 
by  Wehausen  and  Laitone  (1960).  The  differences  occur  as  a result  of 
misprints  in  the  text  of  the  original  reference.  It  is  believad  that 
the  form  of  the  finite  depth  potential  is  correct  a6  given  in  later  sec- 
tions of  this  paper. 

Subsequent  sections  and  appendixes  provide  a formulation  of  the  problem, 
a solution  technique,  a presentation  of  the  effects  resulting  from  pre- 
scribed variations  of  the  problem  parameters  and  a description  and  listing 
of  a computer  program.  The  results  of  the  effort  indicate  that  the  repre- 
sentation and  numerical  evaluation  technique  obtained  herein  will  allow  a 
practical  numerical  calculation  of  the  finite  depth  velocity  potential 
that  can  be  adapted  for  use  with  ship  motion  computer  programs. 
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FORMULATION  OF  THE  PROBLEM 


The  objective  is  to  determine  the  velocity  potential  and  then  the  free- 
surface  elevation  for  the  two-dimensional  problem  of  a pulsating  sourse  of 
arbitrary  strength  fixed  in  a fluid  of  constant  finite  depth.  The  problem 
will  be  formulated  in  terms  of  potential-flow  theory  with  a inear ized 
free-surface  condition.  This  approach  allows  the  following  assumptions: 

1.  The  fluid  is  an  ideal  fluid,  i.e.,  incompressible  and  inviscid. 

2.  Surface  tension  effects  are  negligible. 

3.  The  motion  amplitudes  and  the  fluid  velocities  are 
small  enough  so  that  all  but  the  linear  terms  in 
the  free-surface  and  other  boundary  conditions  may 
be  neglected. 

A complete  discussion  of  linearized  free-surface  theory  is  provided  by 
Stoker  (1957). 

A two-dimensional  coordinate  system  with  the  y-axis  pointing  upward 
and  the  x-axis  in  the  undisturbed  free  surface,  as  shown  in  Figure  1,  is 
used  as  a reference  for  the  problem  formulation  and  numerical  calculations. 


V 


Figure  1 - Reference  Coordinate  System  and  Problem  Geometry 
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where  Q is  the  source  strength,  a is  the  x-location  of  the  source,  b is 
the  y-location  of  the  source,  t is  the  time,  o is  the  pulsating  frequency, 
and  h is  the  depth  of  the  fluid.  The  remaining  terms  in  Equation  [1]  are 
given  by  the  following  expressions: 
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r - /(x-«)2  + (y-b)2 

Tj  ■ /(x-a)2  + (y  4*  2h  + b)2 

is  the  root  of  the  equation:  mQ  tanh  mQh  - v 

The  derivation  of  the  potential  expression  in  Appendix  A follows  a 
method  outlined  by  Thome  (1951)  for  the  infinite  depth  case.  The  form  of 
the  potential  in  Equation  [1]  is  given  by  Uehausen  and  Laitone  (1960) 
with  slight  modification.  The  two  forms  are  equivalent,  however,  as  is 
shown  in  Appendix  A.  The  modification  to  the  Wehauaen  and  Laitone  (1960) 
form  of  the  potential  involves  a change  of  sign  in  one  of  the  terms  within 
the  principal  value  integral.  The  change  is  Indicated  by  an  asterisk  above 
the  modified  sign  in  Equation  [1]. 

The  free-surface  elevation  due  to  a submerged  pulsating  source  was 
selected  as  the  parameter  most  suitable  for  observing  the  influence  of 
fluid  depth,  source  submergence,  and  frequency.  This  was  done  because  the 
free-surface  elevation  is  a more  easily  recognized  physical  parameter  than 
the  value  of  the  potential  itself.  According  to  potential-flow  theory, 
the  free-surface  elevation  is: 

n(x,o;t)  - - j $t  (x,o,t) 

Numerical  determination  of  the  free-surface  elevation  requires  that 
the  potential  be  expressed  in  a form  having  "good"  numerical  behavior,  that 
is,  converging  to  a solution  in  a reasonable  manner  for  all  expected  values 
of  the  physical  parameters. 
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SOLUTION  OF  THE  PROBLEM 


[2] 


Appendix  B provides  a brief  explanation  of  the  process  involved  in 
expressing  the  potential  in  terns  of  a series  expansion  which  can  be  used 
for  numerical  computations.  The  resulting  form  of  the  potential  is: 

*(x’y;t)  “ ^ ~2™h  m°^]  CO>h  My+h>cosh  m^b+lOsin  (mjx-al-at) 

-q)  ■4v2‘H> — 5 — co*  (y+h)cos  m.  (b+h)e-mk^x-*^cos  at 

L,  m.  (hm . +hv  -v) 
k“l  *■  K 

where  m^  are  the  positive,  real  roots  of  the  equation 
tan  ■ -v 

and  mQ  is  the  postlve,  real  root  of  the  equation 
mQ  tanh  m0  ■ v 

The  free-9urface  elevation  thus  becomes: 


n(x,o;t)  - 

“oR 


r(mQ±v 

L 2m0h 


ll 


_moh 


cosh  m, 


h+sinh  2mQh 


cosh  m0(b4h)co8  (mj  x-a|  -at) 


Qa\  "k  + J , . - 

— *— > = t — C09  hcos  m.  (b+h)e 

(hm^  +hv  -v  ) 


'\1  X'£ 


[3] 


sin  at 


A computer  program  has  been  developed  for  numerical  evaluation  of  Equation 
[3].  The  details  of  the  program  are  given  in  Appendix  C. 
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PRESENTATION  AND  DISCUSSION  OF  RESULTS 


The  results  of  a study  that  systematically  varied  the  value  of  the 
principal  parameters,  the  water  depth,  the  depth  of  the  source  below  the 
free  surface,  and  the  frequency  of  oscillation  are  presented  at  the  end 
of  this  section  in  Figures  2-5.  The  computations  for  a given  wave  form 
were  performed  for  constant  values  of  all  the  physical  parameters,  with 
the  exception  of  the  x-distance  of  the  field  point  from  the  source.  The 
x-distance  was  varied  so  that  the  wave  form  in  the  figures  represents  a 
"slice"  of  the  wave  at  an  instant  in  time.  The  source  strength  was  unity 
for  all  computations. 

Figure  2 provides  a comparison  of  the  wave  form  generated  by  using 
the  finite  depth  program  and  that  generated  by  using  a previously  developed 
infinite  depth  program.  The  conditions  for  the  finite  depth  runs  were: 
h (water  depth)  ■ 10  ft 

b (source  submergence)  - -3.0  ft 

o (frequency)  ■ 3.5  rad/sec 

For  the  infinite  depth  runs,  the  source  was  located  at  -3.0  ft  and 
the  frequency  was  3.5  rad /sec.  The  values  of  the  time  parameter  selected 
for  both  programs  were  t - 2.0196,  A. 0392  and  5.3851  sec.  These  values, 
in  combination  with  the  selected  frequency,  result  in  the  ot  term  in 
Equation  [3]  assuming  values  of  ir/4,  tt/2,  and  0 respectively.  Both  phase 
and  amplitude  of  the  resultant  wave  forms  may  be  compared  by  calculating 
the  wave  forms  at  several  different  values  of  time. 
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Figure  2 demonstrate*  that  the  Infinite  and  finite  depth  results  are 
in  close  agreement  for  selected  conditions.  There  is  a small  difference  in 
the  wavelength  (X),  with  the  finite  depth  wavelength  equal  to  16.4996  and  the 
infinite  depth  wavelength  equal  to  16.5158.  Otherwise,  agreement  between  the 
wave  forms  generated  by  the  two  programs  is  quite  close.  Shown  in  Figure  2b 
as  a dashed  line  is  a 5-ft  depth  case  with  other  conditions  the  same  as 
for  the  10-ft  case.  The  results  exhibit  the  expected  characteristics: 
first,  the  wavelength  is  reduced  (X  ■ 15.8938  versus  X ■ 16.4996)  and 
second,  the  amplitude  of  the  wave  is  Increased. 

Figure  3 Illustrates  the  results  of  keeping  water  depth,  frequency,  and 

time  constant  while  varying  the  depth  of  the  source.  For  all  of  the  runs  in 

Figure  3,  water  depth  was  6.0  ft,  frequency  was  3.5  rad/sec  and  time  was 

2.0196  sec  whereas  source  submergence  varied  between  -5.0  and  -1.0  ft.  The 

wavelength  remains  the  same  at  all  values  of  the  source  submergence  because 

wavelength  is  a function  of  m (X  ■ 2ir/m  ),  and  m is  independent  of  the 

o o o 

source  submergence.  The  amplitude  of  the  wave  is  dependent  on  the  source 
submergence,  however,  as  is  the  "range"  of  the  significant  Influence  of  the 
logarithmic  terms.  The  asterisk  in  Figure  3a  indicates  that  the  value  of 
the  term  that  contains  both  the  principal  value  integral  and  the  logarithmic 
terms  is  less  than  0.0001  at  an  x-distance  of  21.0  from  the  source  for  the 
-5.0  source  submergence  case.  The  asterisk  in  Figure  3e  Indicates  that  the 
corresponding  point  for  the  -1.0  source  submergence  case  occurs  at  a some- 
what smaller  value  of  the  x-dlstance  (x  ■ 20.0).  In  addition,  the  importance 
of  the  term  as  a percent  of  the  total  wave  height  is  more  significant  for  the 
b * -5.0  case  than  it  is  for  the  b "-1 .0  case,  thus  the  statement  that  the 
significant  influence  of  the  logarithmic  terms  is  decreased  as  source 
submergence  is  reduced. 
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Figure  4 indicates  the  effect  of  varying  the  water  depth  while  holding 
source  submergence,  frequency,  end  time  constant.  Water  depth  ms  varied 
between  5.0  end  3.0  ft  while  source  subaergence,  frequency,  and  tlae  were 
held  constant  end  equal  to  -2.0  ft,  3.5  rad/sec,  end  2.0196  sec,  respectively. 

The  water  depth  affects  both  wave  aaplltude  end  wavelength.  The 
wavelength  becoaes  significantly  shorter  as  depth  decreases.  In  Figure  4a 
(h  - 5.0)  the  wavelength  is  15.8938  whereas  in  Figure  4c  (h  - 3.0)  the  wave- 
length is  14.3058.  The  wave  aaplltude  Increases  as  water  depth  decreases. 
Also,  as  water  depth  tends  toward  small  values,  the  rate  of  increase  in  wave 
height  becoaes  large. 

Figure  5 shows  the  wave  fora  that  results  froa  a lower  frequency  of 
oscillation.  Part  a of  Figure  5 corresponds  to  part  a of  Figure  4 in  all 
paraaeters  except  frequency.  The  wave  fora  in  Figure  5a  is  for  a frequency 
of  2.5  whereas  the  wave  fora  in  Figure  4a  is  for  a frequency  of  3.5.  The 
saae  correspondence  exists  between  Figures  4b  and  5b. 

The  lower  frequency  results  in  longer  wavelengths  and  a slightly 
saaller  aaplltude  coapared  with  the  higher  frequency  case.  As  in  previous 
examples,  the  results  shown  in  Figure  5 Indicate  that  as  water  depth 
decreases,  wavelength  decreases  and  wave  height  Increases.  For  the  same 
change  in  paraaeters,  however,  the  effects  of  depth  on  wavelength  are 
slightly  more  noticeable  for  the  lower  frequency  than  for  the  higher 
frequency  results.  For  example,  the  ratio  of  the  wavelengths  in  Figures 
4a  and  4b  was  1.04  whereas  the  corresponding  ratio  for  Figures  5a  and  5b 
was  1.08.  A comparison  of  the  amplitudes  in  Figures  4s  and  5a  shows  that 
the  wave  aaplltude  for  the  lower  frequency  was  less  than  that  for  the  higher 
frequency  (0.0472  versus  0.0460). 
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CONCLUSIONS 


As  a result  of  this  study,  the  following  conclusions  may  be  drawn: 

1.  The  techniques  employed  In  this  study,  for  a numerical  evaluation 
of  the  finite  depth  potential,  will  allow  the  extension  of  the  source 
distribution  method  to  provide  two-dimensional  added-mass  and  damping 
characteristics  in  a fluid  of  finite  depth  for  use  with  modern  ship  motion 
strip  theories. 

2.  Comparisons  between  wave  forms  computed  from  both  finite  and 
Infinite  depth  potentials  indicate  that  fluid  depth  does  not  begin  to 
significantly  alter  the  wave  form  until  the  ratio  of  wavelength  to  water 
depth  is  1.7  or  greater. 

3.  The  deeper  the  submergence  of  the  source,  the  lower  the  amplitude 
of  the  wave  form.  Wavelength  is  not  affected  by  source  submergence.  Also, 
as  the  source  is  placed  closer  to  the  free  surface,  the  significant  effects 
of  the  logarithmic  terms  in  the  potential  expression  decay  more  rapidly 
with  distance  from  the  source, 

4.  As  fluid  depth  decreases,  wave  amplitude  increases  and  wavelength 
decreases.  The  rate  of  Increase  in  wave  amplitude  is  increased  as  fluid 
depth  decreases. 

5.  Pulsating  frequency,  or  frequency  of  oscillation,  affects  both 
wavelength  and  wave  amplitude.  With  other  parameters  held  constant,  the 
effect  of  fluid  depth  on  wavelength  becomes  more  noticeable  and  wave  amplitude 
is  reduced  as  frequency  is  decreased. 
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APPENDIX  A 


i 


DERIVATION  OF  THE  VELOCITY  POTENTIAL  FOR  WATER  OF  FINITE  DEPTH 

The  derivation  of  the  velocity  potential  for  a aource  puleating  in  a 
fluid  with  free  surface  and  finite  depth  le  presented  in  this  appendix. 

The  derivation  follows  a method  suggested  by  Thorne  (1953)  for  the  infinite 
depth  case.  However,  the  computer  program  used  to  evaluate  the  potential 
is  based  on  a form  due  to  John  (1950)  and  elaborated  on  by  Uehausen  and 
Laitone  (1960).  The  two  forma  are  equivalent,  however,  as  will  be  shown 
at  the  end  of  this  appendix. 


COORDINATE  SYSTEM  AND  IDENTIFICATION  OF  PHYSICAL  PARAMETERS 

Figure  A.l  shows  the  coordinate  system  used  for  the  derivation  and 
illustrates  the  basic  geometric  parameters.  The  y-axls  is  positive 


downward  and  the  x-axis  is  in  the  undisturbed  free  surface  of  the  fluid. 


Vb)  ry 


Free  Surface 


' • O'.* ) 


Bottom 


Figure  A.l  - Coordinate  System  and  Problem  Geometry  for  Derivation 


of  the  Potential  Expression 


i 

I 
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The  coordinate  system  is  the  same  as  employed  by  Thorne  (1953) 
but  different  than  that  employed  by  Wehausen  and  Laitone  (1960).  The 
latter  employ  a system  that  has  the  y-axis  positive  in  the  upward  direction, 
as  shown  in  Figure  1.  The  computer  program  also  employs  the  coordinate 
system  shown  in  Figure  1 and  not  the  one  given  in  this  appendix. 

The  physical  parameters  are  define^  as  follows: 
h - depth  of  fluid  as  measured  from  the  mean  free  surface 
to  the  bottom;  the  bottom  is  assumed  to  be  uniform 
a - x-location  of  the  source 
b - y-location  of  the  source 
0 - tan  1 (x-a)/(y-b) 

0'  - tan  1 (x-a)/(y+b) 
r - /(x-a)^  + (y-b)^ 

r ' ■ »4x-a)^  + (y+b)^' 

x - x-location  of  the  point  at  which  we  wish  to  compute 
the  value  of  the  potential 

y - y-location  of  the  point  at  which  we  wish  to  compute 
the  value  of  the  potential 
o - frequency  of  oscillation  of  the  source 

STATEMENT  OF  PROBLEM 

Since  this  analysis  is  for  the  two-dimensional  case,  the  singularity 
will  be  logarithmic  in  form.  The  statement  of  the  problem  has  been  presented 
in  the  main  text,  but  it  is  repeated  here  for  continuity  of  presentation. 
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and  because  of  the  possible  confusing  over  the  different  coordinate  systems. 
We  seek  a function  of  the  form 

4(x,y;t)  “ 4>1(x,y)cos  at  + <|>2(x,y)sin  at  [ A— 1 ] 

defined  for  y >,  o,  except  at  (a,b)  and  satisfying  the  conditions 

1.  V^<p ^ - 0,  i - 1,2,  except  at  (a,b) 

^i  2 

2.  + — - 0,  i - 1,2,  v - a /g 
<*4., 

3.  -gy  - 0,  i - 1,2,  at  y - h 

4.  The  radiation  condition,  i.e.,  the  waves  are  outgoing  as  x tends 
to  infinity 

5.  4j_(x,y)  - log  r/r'  + <J>q (x ,y) 


DERIVATION  OF  POTENTIAL 


Thorne  (1953)  suggested  the  following  form  for  the  <J>q (x ,y ) potential: 


*0(x,y)  - J [ago  sinh  ky+B(k)cosh  k(h-y)Jcos  k(x-a)  dk 

This  is  a harmonic  function  that  satisfies  condition  1 since 

2[ 


[A-2] 


oxx 


oyy 


■h 


A(k)sinh  ky+B(k)cosh  k(h-y) 
A(k)8inh  ky+B(k)cosh  k(h-y) 


cos  k(x-a)dk 
cos  k(x-a)dk 


and  hence  4 +4  -0.  The  4_(x,y)  potential  will  be  determined  later 

oxx  oyy  2 

to  satisfy  the  radiation  condition  (Item  4) . The  function  A(k)  and  B(k) 
must  be  chosen  such  that  the  integral  in  [A-2]  converges  in  the  region 
-oo  < x < oo , 0 <_  y h 

To  accomplish  this  task,  the  expression  for  ^(x.y)  is  put  into  the 
free-surface  condition  (Item  2)  to  give: 
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m 

—■  log  + vlog  |k(k)slnh  ky+B(k)eoah  k(h-y^cos  k(x-a)dk 

+ vJ|A(k)*lnh  ky+B(k)coeh  k(h-y)Jcos  k(x-a)dk  ■ 0 

At  y ■ o,  the  following  la  true  from  conelderatlon  of  the  geometry  of 


t A—  3 1 


the  problem: 


3 . r -2b 

— log  - 5 2 

3y  r (x-a)Z+bZ 


(A-4  ] 


The  derivative  In  [A— 3]  may  be  expreased  aa  followa 


oo 

Hi 

•It 


A(k)slnh  ky+B(k)cosh  k(h-y)J  cos  k(x-a)dk 
kA(k)cosh  ky-kB(k)elnh  k(h-y)Jcos  k(x-a)dk 


[ A— 5 ] 


Substituting  the  fA-4]  and  [ A— 5 ] relations  into  [A— 3]  for  y ■ o gives 


lkA(k)-B(k) [k  alnh  kh-vcosh  kh]|cos  k(x-a)dk  • 


2 2 

(x-a)N-bZ 


[ A-6  J 


The  "well-known"  Identity* 


r — g 

I e cos  royrlx  - — ^ — =-  for  B > 0 

1 

can  be  written  in  terms  of  the  problem  variables  as 

2 fe  ^cos  k(x-a)dk  • —s — r 

J (x-a)z+bz 

o 

Using  this  identity  with  Equation  [A-6]  and  equating  the  Integrands  of  the 
resulting  Integral  expressions  gives  the  following  equation  in  A(k)  and  B(k): 


* See  Hodgmen  (1959),  page  314,  number  430. 
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BOO  - -fcAQ0+2«~kb 
1 ' -kslnh  kh+vcosh  kh 


[A-7] 


Consider  now  the  bottom  condition  (Item  3).  The  final  form  of  the 
function  must  satisfy  this  condition  for  y ■ h.  At  y • h,  the 
problem  geometry  gives : 


3 . r 

W 108  r* 


h-b 


(x-a)2+(h-b)2  (x-a)T+(h+bj2 


h+k 


•ft 


-k (h-b)  -k (h+b) 
“6 


cos  k(x-a)dk 


[A-8] 


Setting  y ■ h in  Equation  [A-5]  and  equating  the  right-hand  side  of  this 
expression  to  the  Integral  in  [A-8]  yields  the  integral  equation 


1- 


kA(k)cosh  khcos  k(x-a)dk 


■4 


e khcosh  kbcos  k(x-a)dk 


[A-9  ] 


It  can  now  be  observed  chat  in  order  to  satisfy  the  bottom  condition  the 
following  must  be  true: 

-kh 


A(k)  - 


-2e 


kcosh  kh 


sinh  kb 


[ A-10 ] 


The  combination  of  Equations  [A-7]  and  [A-10]  provides  an  explicit 
expression  for  B(k): 


nfifl  - 2(e  ktlslnh  kb-f-e  kbcosh  kh) 
cosh  kh  (kslnh  kh-fvcosh  kh) 


[ A— 1 1 ] 


Introducing  the  expressions  for  A(k) , [A-10],  and  B(k) , [A-ll],  into 
Equation  [A-2]  gives  the  desired  expression  for  the  potential: 


*„(x,y) 


cosh  k(h-b)cosh  k(h-y) 
cosh  kh(vcosh  kh-ksinh  kh) 


-kh 


e sinh  kbslnh  ky 
kcosh  kh 


cos  k(x-a)dk 

[ A— 12  ] 
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Th«  next  step  la  to  detaralna  d2(x,y),  whlch  *■  tha  aacond  tarn  in  tha 
potantlal  axpraaalon  In  [ A— 1 ] ; $2(x,y)  ouat  ba  aalactad  auch  that  +(x,y;t) 
repraaanta  an  outgoing  vava  aa  x tanda  to  Infinity  (tha  radiation  condition). 
To  accompllah  thla  atap,  tha  form  of  the  potantlal  for  larga  valuaa  of 
x must  be  found.  Thla  requlrea  avaluatlon  of  the  Integral  In  [A-12]  around 
a contour  formed  by  the  quadrant  R(k)  > 0,  I(k)  > 0 but  excluding  any  singu- 
larities  that  occur  along  the  R(k)  axis.  Figure  A. 2 shows  such  a contour. 


R(k) 


Figure  A. 2 - Contour  of  Integration 
It  is  convenient  to  rewrite  the  integral  express  for  $Q(x,y)  differently 
than  that  shown  in  [A-12]  to  give: 


00 

fi 


— k.h 

kcosh  k(h-b)cosh  k(h-y)-e  slnh  kbslnh  ky(vcoah  kh-kalnh  kh) 
kcosh  kh  (vco*h  'kh-ksihh 'kh) 


cos  k(x-a)dk 


.cos  k(x-a)dk 


[ A— 13] 


We  are  concerned  with  the  case  where  k Is  greater  than  zero.  The  only 


singularity  in  the  integral  expression,  [A-13],  therefore,  occurs  at  tha 
value  of  k - mo,  where  mQ  la  the  poatlva,  real  root  of  the  equation: 


vcosh  kh-kslnh  kh  ■ 0 


[A-14] 
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The  integral  expression  in  [A-13]  may  now  be  written  as  the  sum  of  the 
integrals  along  the  parts  of  the  contour  of  integration  plus  the  value 
of  the  residue  term  at  k ■ n0: 


where  along  the  axes 


'18Xd8 

^ ^ 

-iff [ Res idue ] eim°x 


■jj$s  ■•nXd" 

+ fltSsifi.  a11>«10-d(Re19)] 
q(Reiy)  J 

o 


i 

o 


[A-15 


k - 6+iri , dk  ■ d6  for  i 

dk  • idn  for  J 


and  along  the  arc 


k - Re10,  dk  - iRei0dO 


R > 0,  0 <_  0 <_  ff/2 


The  individual  expressions  in  [A-15]  are  now  examined. 
The  integral  along  the  arc 

* 4?  IQ  10 

J ,(«.le) 
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tends  toward  zero  as  R -*•  *. 

The  integral  along  the  imaginary  axis 


a 


0 

tends  toward  zero  as  x -*•  ». 

The  integral  along  the  real  axis  is  the  principal  value  integral  we  seek  to 
evaluate.  The  integral  on  the  left-hand  side  of  [A-15]  is  an  integral  about 
a closed  analytic  region  and  is  zero  by  the  Cauchy  integral  theorra.  The 
value  of  the  residue  term  is  determined  in  the  following  manner: 

Define 


f(k)  ■ n'.V  and  b.  - the  residue  at  k - m 
q(k;  1 o 

According  to  Churchill  (1960) , the  residue  of  f(k)  due  to  a singularity  at 
k - mo  exists  if 

a.  p (mo)  4 0 

b.  q (mo)  4 0 

c.  q' (mo)  4 0 

Since  mQ  is  defined  as  the  root  of  the  q(k)  function  and  q(mo)  ■ 0, 
condition  b is  satisfied.  Likewise  condition  a is  satisfied  since  p(m0)  4 
0 if  k > 0,  and  such  a restriction  on  k has  already  been  specified.  The 
expression  for  q'(k)  must  be  examined  to  determine  whether  condition  c is 
met  by  the  f(k)  function.  Now,  q'(k)  may  be  expressed  as: 

2 

q’(k)  » vcosh  kh  + 2kvhsinh  khcosh  kh 
- ksinh  2kh  - k^hcosh  2kh 
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Alio,  at  k ■ m0, 

v " "V>  co'ah"m^h  and  q*  (m<,)  - ^ >inh  n<)h  ♦ a<J2h 
In  order  to  aatlafy  condition  c the  expression 


slnh  2ra0h  -f  2m0h  ^ 0 


lA-16] 


must  be  true.  The  only  root  of  Equation  [A-16]  occures  at  2moh  - 0 
or  at  mo»0  since  h is  always  greater  than  zero.  The  value  of  m,, 
cannot  be  zero  If  h is  greater  than  zero  and  therefore  condition  c is 
satisfied.  The  residue  term  in  Equation  [A-15]  thus  becomes 


, m -2cosh  mfi(h-b)cosh  ny, (h-y) 

<1 L * _ J —l_  ^ L 


2m0h  + slnh  2m0h 


[ A- 1 7 ] 


All  of  the  integrals  and  terms  In  Equation  [A-15]  have  now  been  evaluated. 

Taking  the  real  part  of  the  resultant  form  of  the  integral  yields  an 
expression  that  is  Introduced  Into  the  potential  function,  Equation  [ A— 12 ] , 
to  yield  the  following  form  for  $ (x,y)  for  large  x: 


* (x,y)  - i’cf,h ■taft-Ms;*)! -SoC-y?  ,1„  v 
o 2mch  + slnh  2moh  ° 


2moh 


[A- 18] 


As  x becomes  large,  the  logarithmic  terms  dissipate  rapidly.  This  means 
that  the  stun  of  $o(x,y)  and  $2(x*y)  must  combine  to  ensure  that  the  radiation 
condition  (Item  4)  is  satisfied.  A function  that  performs  this  task  Is: 


-ATTcosh  mf,(h-b)cosh  ty, (h-y) 
2moh  + slnh  2moh 


cos  mQx 


[ A— 19 ] 


All  of  the  functions  required  to  obtain  the  total  velocity  potential 
expression  are  now  available. 

Introducing  ♦ (x,y),  Equation  [A-12],  $_(x,y),  Equation  [A-19], 

O c 

and  the  logarithmic  terms  into  Equation  [A-l]  provides  an  expression  for 
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the  velocity  potential  that  satisfies  all  requirements  and  constraints: 


$(x,y;t) 


00 

JL  1 L + ijr  [cosh  k(h-b)cosh  k(h-y) 

2tt  r'  j | cosh  kh(vcosh  kh-ksinh  kh) 

° 

e sinh  khsinh  ky]  , , s „ 

: r — — ^Jcos  k(x-a)dk  cos  at 

kcosh  kh 

_ 2qcosh_qn  (hjd? ) cqsh  _ma(_h-y )_  coa  (x..a)atn  ot 
2m0h  + sinh  2m0h  ° 


[A-20] 


This  is  the  same  as  the  form  derived  by  Thorne  (1953) . 


ALTERNATE  FORM  FOR  THE  VELOCITY  POTENTIAL 

The  most  convenient  form  for  the  velocity  potential  is  not  the  one 
expressed  by  Equation  [A-20].  Because  of  planned  adaptations  of  the 
numerical  evaluation  of  the  potential,  a form  given  by  Wehausen  and  Laitone 
(i960)  is  more  desirable.  The  principal  differences  between  the  two 
alternate  forms  stem  from  the  reference  coordinate  systems  employed  by  the 
individual  investigators.  As  mentioned  previously,  Thorne  (1953)  employs 
a coordinate  system  with  the  y-axis  positive  in  the  downward  direction,  as 
shown  in  Figure  A.l.  Wehausen  and  Laitone  (1961),  however,  base  their 
potential  on  a coordinate  system  with  the  y-axis  positive  in  the  up 
direction  as  shown  in  Figure  1.  Other  differences  are  due  to  slightly 
differing  methods  of  representing  geometric  parameters.  To  convert  from 
the  potential  form  in  Equation  [A-20]  to  the  equivalent  expression  as 
given  by  Wehausen  and  Laitone  (1960)  requires  the  following  transformations 
and  identities: 
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y - -y,  x - x 

a - a , b = -b 


h - h 


Ct 

] 

J 


i , 2jLa2,l/2 

1 .--e.SS-J£  e'8udy  - log  , 6 > 0 

u 8 


-Bu  -hy 


e -e 


dy  - log  - , 8 > 0 


e m°hsinh  nnh  2e  ^^coah  ny,h  m mp-v 
vh+.lnh^h  ' 2"°h-‘lnh  2moh  h«/-hv2+v 


[A-21] 


Substitution  of  [A-21]  into  [A-20]  and  algebraic  manipulation  results  in 
the  desired  form  of  the  velocity  potential: 


Of,  r . , r2  _ I I k+v  e khcosh  k(h+b)cosh  k(y+h) 

$(x,y;t)  2n  [_  °8  h h ] ( k ksinh  kh-vcosh  kh 

-kh) 

+ V~}  dk} 


cos  at 


cos  k(x-a) 

[ A— 22 ] 


y+m n e m°hsinh  mnhcosh  mr,(h4’b)cosh  mn(y+h)  cos  ny,(x-a)  g^n 
ra°  vh  + sinh2  mph 


/ 9 7 

where  r - /(x-a)  + (y-b)* 

/ 2 ? 

r2-  /(x-a)  + (y+2h+b)  , and  y is  positive  upward. 

The  source  is  thus  located  at  a distance  y-b  beneath  the  free  surface  and 
the  bottom  is  located  at  a y«-h  distance  from  the  free  surface. 
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APPENDIX  B 


SERIES  EXPANSIONS  USED  TO  NUMERICALLY  EVALUATE  THE  VELOCITY  POTENTIAL 


A brief  outline  of  the  procedure  used  to  derive  a series  expansion  for 
the  finite  depth  velocity  potential  is  provided  in  this  appendix.  A more 
detailed  explanation  may  be  found  in  John  (1950). 

The  integral  expression  for  the  potential  is  given  by  equation  [ A— 2 2 ] 
but  is  repeated  here  for  convenience: 


*(x,y ;t) 


2ir 


-kh 


log  - + log 


— _ir  I k'fv  e " “cosh  k(h+b)cosh  k(y+h)cos  k(x-a) 
h j I k ksinh  kh  - vcosh  kh 


-kh 

-r  > dk 


cos  at 


[A— 22 ] 


qVfmo  e m°^slnn  mphcosh  mp (h+b) cosh  mo(y+h)cos  mn(x-a)  sinot 
m°  vh  + sinh^  mQh 

The  first  step  is  to  express  the  coefficient  of  the  sin  at  term  in  the 
potential  expression,  equation  [A-22],  in  a different  form  through  use  of 
the  following  identity: 


e m°^slnh  mph  _ mp-v 

2 *22 
vh  + sinh  mph  hirip  -hv  +v 

Next,  the  identities  provided  by  equations  [A-21]  are  used  to  include  the 
logarithmic  terms  under  the  integral  in  the  coefficient  of  the  cos  at  term 
in  equation  [A-22].  The  integrand  of  the  resulting  integral  is  now 
expressed  in  the  form  of  a partial  fraction  expansion  and  integrated  term 
by  term.  The  resultant  expression  for  the  potential  is: 
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[B-l] 


$(x,y;t)  - -fi-  -~-v — r—  cosh  m (y+h)coah  m (b+h)ain  (a  |x-a|-ot) 

m_  , L . L o o o 


o hm0  -hv 


2 2 
+ w 


- Q j cos  mk(y+-h)cos  mk(b+h)e-mk^  x-aUos  ot 

k-1  \ hmk  + hv  “v 


where  mQ  Is  the  positive  real  root  of  the  equation 

m tanh  m h - v 
o o 

and  are  the  positive,  real  roots  of  the  equation 
tan  m^h  - -v 


The  asterisk  under  the  cos  ot  In  equation  [B-l]  indicates  that  this 
term  differs  from  the  corresponding  term  in  the  expression  given  by 
Wehausen  and  Laltone. 


An  expression  for  the  free-surface  elevation  may  now  be  found  by 
recalling  that 

n(x,o  ;t)  “ ” ~ 4>t  [B-2] 

Introducing  the  time  derivative  of  the  potential  expression  in  equation 
[B-2]  results  in  the  following  equation  for  the  free-surface  elevation: 


n(x,o;t) 


2 2 

_ 171  -V 

Qo  o 


gm0  hmQ  -hv  +v 


— cosh  m^hcosh  mo(b+h)cos  (mo|x-a|-ot) 


[B-3] 


- 9s. 

g 


Ei 

— --  7-7-  cos  rn^cos  m^b+h) 


k-1 


2 2 
+v 

1 

2 L 

hm^  +hv  -v 


e-“k> 


x-a 


sin  ot 
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A slightly  different  form  of  this  expression  was  adopted  for  use  with  the 
computer  program.  This  form  employs  the  identity 


2e  ^ cosh  nv,h 


2m0h+sinh  2moh  hm^-hv^v 


to  change  the  form  of  the  first  term  in  equation  [ B— 3 ] . The  final 
expression  for  the  free-surface  elevation,  as  programmed,  is: 


,(x  0.t)  . 1 32.  On0-Hv)e~m°h  cosh2 
' m0g  2m0h+sinh  2moh 


cosh  mQ(b+h)cos  (m0|x-a)  -at) 


\ ’ m,  +v  -m, 

~ i)  ^ v,  2 .1 C0S  VC0S  %(b+h)e 

,1c  hm,  +hv  -v 

k-1  k 

This  representation  is  valid  for  b>0,  h>0,  and  r + 0. 


-n^lx- 


sin  at 


i 


APPENDIX  C 


COMPUTER  PROGRAM,  LISTING,  SAMPLE  INPUT  AND  SAMPLE  OUTPUT 

The  computer  program  has  been  assembled  for  use  with  a CDC  6700 
computer  and  is  written  in  the  Fortran  IV  language.  The  program  is  easily 
adaptable  to  any  modern  high-speed  computer  with  only  minor  changes  in  the 
Fortran  IV  program  listing.  No  special  input,  output,  or  external  storage 
devices  are  required.  The  program  consists  of  one  main  program  and  two 
subroutines.  The  total  program  is  entitled  FINDEP;  a listing  is  provided 
at  the  end  of  this  appendix. 

MAIN  PROGRAM 

The  main  program  performs  all  input,  output,  and  most  of  the  numerical 
calculations  involved  in  determining  the  free-surface  elevation.  The 
potential  itself  is  not  actually  determined,  but  the  program  could  be 
easily  modified  to  provide  such  data. 

The  infinite  sum  term  found  in  equation  [B-4]  is  evaluated  to  N number 
of  terms  in  the  main  program.  A convergence  criterion,  stated  below, 
determines  the  value  of  N for  each  value  of  x at  which  computations  are 
performed.  For  computational  purposes,  the  value  of  all  parameters  were 
considered  constant  except  the  x-distance  of  the  field  point  from  the  source. 
The  resultant  wave  form,  therefore,  represents  a "slice"  of  the  wave  at  an 
instant  in  time. 

The  convergence  criterion  of  the  infinite  series  is  baaed  on  the 
difference  between  the  average  values  of  the  summation  term.  Figure  C.l 
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shows  a typical  plot  of  the  term  PHl2(k)  versus  k for  x * a - 0.  The 
symbol  PHI2(k)  represents  the  sum  of  the  individual  terms  in  the  infinite 
series.  As  can  be  observed  in  the  figure,  convergence  of  the  series  to  a 
representative  value  is  quite  slow.  An  estimate  of  the  convergence  value 
is  made  by  determining  the  maximum  and  minimum  values  of  PHI2(k),  taking 
the  average  of  these  numbers,  and  returning  the  value  of  the  last  average 
when  the  difference  between  previous  averages  is  less  than  or  equal  to 
0.0001.  The  value  of  N at  which  the  criterion  is  met  is  the  number  of  terms 
at  which  the  series  is  truncated.  If  the  convergence  criterion  is  not  met 
by  the  time  N ■ 189,  the  average  of  the  last  minimum  and  last  maximum  value 
is  returned.  Thus,  the  value  of  N never  exceeds  189. 

For  values  of  (x-a)  not  equal  to  zero,  the  series  convergence  is 
greatly  accelerated  because  of  the  rapid  decay  of  the  exponential  term. 
SUBROUTINE  R00T1 

The  R00T1  subroutine  computes  the  root,  n,  of  the  transcendental 
equation 

n tanh  n ■ B 

where  n - m^h  and  B ■ vh  The  technique  used  to  solve  the  equation  employs 
the  "Newton-Raphson"  iteration  technique  with  a convergence  criterion  that 
requires  the  Incremental  estimate  for  the  root  to  be  less  than  or  equal  to 
0.C0001.  The  "Newton-Raphson"  technique  works  quite  well  and  a suitably 
accurate  root  is  generally  attained  within  seven  iterations.  Figure  C.2 
presents  a representative  graphical  solution  to  the  equation  for 
hypothetical  problem  parameters.  It  is  clear  from  the  figure  that  as  B 
gets  large,  the  value  of  B/n  approaches  1.0,  and  that  the  root  will 
always  be  positive  and  real.  A detailed  description  of  the  "Newton-Raphson" 

technique  is  given  in  Hildebrand  (1956) . 
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SUBROUTINE  ROOT2 


The  ROOT2  subroutine  computes  N number  of  roots  of  the  transcendental 
equation 

s tan  8 • -B 

where  s ■ m^h  and  B ■ vh.  As  with  subroutine  R00T1,  the  "Newton-Raphson" 
Iteration  technique  is  employed  to  determine  the  value  of  each  of  the 
required  N number  of  roots. 

Figure  C.3  provides  a representative  plot  of  a graphical  solution  to 
the  subject  equation  for  typical  values  of  the  problem  parameters.  The 
positive,  real  roots  of  the  equation  are  sought,  and  the  roots  must  there- 
fore lie-  in  the  second  and  fourth  quadrants  where  the  tangent  function  is 
negative.  As  s gets  large,  the  value  of  the  root  approaches  kir,  where  k 
is  the  summation  index  as  shown  in  equation  [B-4]. 

INPUT 

The  required  input  is  identified  below  together  with  the  data  card 
requirements.  Two  data  cards  are  required. 

Data  Card  1 (7F10.4) 

The  following  floating  point  parameter  values  provide  input  via  data 
card  1.  They  are  listed  in  the  order  in  which  they  must  appear  on  the  card. 


1. 

Q - 

source  strength 

4. 

H - 

water  depth 

2. 

G - 

gravitational  constant 

5. 

A ■ 

x-location  of  source 

3. 

F - 

frequency 

6. 

B - 

y-location  of  source 

7. 

T - 

time 
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Data  Card  2 (F10.4,  16) 

The  parameters  that  provide  Input  via  data  card  2 are  as  follows: 

1.  STEPX  ■ Incremental  value  of  x 

2.  STEPY  ■ Incremental  value  of  y 

3.  STEPZ  - number  of  points  at  which  the  free-surface 

elevation  will  be  calculated  plus  1 

The  size  of  the  increment  step  within  the  range  of  x is  determined  by 
STEPX.  The  value  of  STEPY  must  always  be  zero.  The  value  of  (STOPZ-1) 
determines  the  range  of  the  calculation.  Thus,  the  number  of  values  of  x 
at  which  the  elevation  will  be  calculated,  M,  is 

(STOPZ-1) *STEPX  - M 

For  the  present  intended  use  of  the  program,  inclusion  of  an  iteration 
on  y is  not  necessary.  It  was  Included,  however,  because  of  planned 
extensions  to  the  program  that  include  the  calculation  of  the  velocity 
potential  at  points  other  than  on  the  free  surface. 

The  program  is  constructed  such  that  all  calculations  start  at  a 
value  of  x - 0.0. 

SAMPLE  INPUT  AND  OUTPUT 

A sample  set  of  parameter  values  and  constants  is  identified  below  in 
terms  of  the  variables  in  equation  [B-4] . 

Data  Card  1 Q - 1.0,  g - 32.2,  b - 3.5 

h - 6.0,  a - 0.0,  b - -4.0,  t - 2.0196 

Data  Card  2 Incremental  value  of  x ■ 1.0 


Incremental  value  of  y “ 0.0 
Range  of  x values  +1-42 


The  results  of  operating  the  program  with  the  above  input  values 
are  shown  in  the  "OUTPUT"  listing  at  the  end  of  this  appendix.  Other 

Chan  titles  and  ldentif ications , the  first  information  output  is  a listing 
of  the  input  data  Certain  interim  calculation  terms  are  printed  next: 


SIGMA 

• 

2. 

c / g 

GMU 

m 

mo 

ANGL1 

m 

cosh  m0(b+h) 

SINF 

m 

sin  ot 

C0EF1A 

m 

m^+v  -mnh 
— “ — e ° cosh  nuh 
mo  ^ 

C0EF1B 

C0EF1 

cosh  mr,(b+h) 

2moh  + sinh  2m0h 
2*C0EFF1A*C0EF1B 

and  the  calculated  wavelength.  The  final  output  list  includes  the  listing 
of  the  wave  height  as  a function  of  x and  y,  the  number  of  loops  plus  1 
(N  + 1)  required  for  the  infinite  series  to  meet  the  convergence  criteria 
(LOOPS),  the  last  value  of  m^  for  each  x iteration  (M(k)),  an  interim 
coefficient  value  (C0EF2),  the  coefficients  of  the  cos  (n^l  x-a j -ot)  term 
in  Equation  [2]  (PH11),  and  the  coefficient  of  the  sin  ot  term  in  Equation 
[B-4]  (PHI2A(k)) . 

The  execution  time  for  the  program  is  in  the  neighborhood  of  10  sec 
for  the  sample  problem.  The  greatest  amount  of  time  is  consumed  in 
determining  the  convergence  value  of  PHl2(k)  for  the  x - a ■ 0 case.  The 
program  is  quite  economical  to  operate. 


PROGRAM  LISTING 


PROGRAM  FINDER < INPUT • OUTPUT » T APE5« INPUT t T APE6-0UTPUT > 

DIMENSION  PHI? (200) tX (200) ,Y(?©0>  »PHll (200) »ETA(?00) ♦N0(?00) » APG (? 
100) « ARGH ( 200 ) «PMU< 200) tBMU(200) *WAYOUT(200> «COEF2(?00) *PMIMA(100) , 
?PHI?A(100) 

REAL  LAMBDA 
INTEGER  STOPZ 

COSH (W) »0«5# (EXP (W) »1 • 0/EXP <W) ) 

SINH(N)»0.5*(EXP(N)-1.0/EXP(W) ) 

READ (5 ♦ 1 ) Q.G.F.H.A.BtT 

1 FORMAT (7F10»4) 

RE AO (5*2)  STEPX*STEPY,STOP? 

2 FORMAT (?F10*4«I6) 

SIGMA«F«F/G 

SIGMA?sSIGMA*SIGMA 

SIGMAH»SIGMA*H 

CALL  ROOT1 (5IGMAH»FMU»N) 

GMU>FMU/H 

LAMBOA»6.2B32/GMU 

GMU2=GMU*GMU 

FRFQ*F#T 

POWA»EXP(-FMU) 

AMPLI=0*F/G 

COSF«COS(FREO) 

SINF*SIN(FREO) 

ANGL1«C0SH(GMU*(B^H> ) 

ANGL2*SINH(2.«FMU) 

ANGL3*C0SH(FMU) 

COEF 1 A* ( GMU*S I GMA ) *PONA* ANGL  3/GMl  I 
COEF 1 B* ANGL 1 / ( 2 . *FMU ♦ ANGL2 ) 

COEF  1 «2  . *COEF  1 Ar*COEF  1 B 

X ( 1 ) *-STEPX 

Y(1)»-STEPY 

DO  100  JJ«2tSTOPZ 

X ( JJ) »X ( JJ-1 ) ♦STEPX 

Y(JJ)*Y(Jj-l)*STEPY 

ARG ( JJ) «GMU* ( ABS ( X ( JJ) -A) ) -FRFQ 

ARGH ( J J) aGMU* (Y ( JJ) *H) 

PHI 1 ( J J ) *COEF 1 *C0SH ( ARGH ( JJ ) ) 

PHI2(1)«0. 

COEF2 ( 1 ) »0. 

KK*0 

KMIN*0 

KMAXaO 

DO  90  K»2* 190 
L«K 

CALL  R00T2 (SIGMAH* ANU* I *K) 

IF (I. GT. 100)  GOTO  50 
BMU  ( K ) * AMU 
AMU«AMU/H 
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APGO-AMU*(Y(JJ)«H) 

ARG1*AMU*(R*H) 

POW- AHU«ABS ( * ( JJ) -A ) 

AMU?«AMU*AMU 

C0EF2 (K) • < AMU?*SIGMA2) / ( AHt« <M«AMU/>*H*SIGMA?-SIGMA) ) 

If (POW.GE. 200.0)  GO  TO  56 

PM ! 2 <K ) *PMI  2 <K- 1 ) ♦COEF2 ( K ) «COS  < APG 1 ) *COS ( ARGO ) *EXP  f -POW) 

GO  TO  57 

56  PMIN*PHI2(K-1) 

PMAX*PHI2(K-1) 

GO  TO  A3 

57  IF(K.GT.4)  GO  TO  40 
GO  TO  90 

40  IF(PHI2(K-2).E0.PMI2<K-1).AND.PHI2(K-1).E0.PMI2(K))  GO  TO  39 
IF<PHI2(K-?) .GE.PMI2 (K-l) . AND.PHI2 (K-l ) .LE.PMI20O ) GO  TO  41 
lFCPHI2(K-?).LE.PMl2(K-l).AN0.PHI2(K-l).GE.PMI2(Kn  GO  TO  42 
GO  TO  90 

41  KMIN«KMIN#1 
PMIN«PHI2(K-1) 

IF(KMAX.EO.O)  GO  TO  90 
GO  TO  43 

42  KMAX«KMAX«1 
PMAX«PH!2(K-1) 

IF(KMIN.EO.O)  GO  TO  90 
GO  TO  43 

39  PMIN«PMI2(K-1) 

PMAX«PHI2 (K-l ) 

43  KK»KKM 

PMIMA(KK)*.S*(PMIN*PMAX) 


22  IF<KK.GE.2>  GO  TO  44 
GO  TO  90 

44  OIFF*PMIMA(KK)-PMIMA(KK-l) 

IF(ABSCDIFF) .LE..0001)  GO  TO  95 
90  CONTINUE 
95  NO ( JJ) «L 

PMI2(N0(JJ))*PMIMA(KK) 

ETA ( JJ) *AMPL I* (PHI 1 ( JJ) *COS ( ARG < JJ)  ) “PH 1 2 (NO  I JJ) ) *SINF ) 

PMU(JJ)*AMU 

WAYOUTCJJ)»COEF2CL) 

PHI2A(JJ)*PMI«A(KK) 
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100  CONTINUE 
GO  TO  30 
50  WRITE (6.55)  I 
55  FORMAT(20X.3H  I*. 16) 

WRITE (6*58)  JJ.K 

58  FORMAT (20X.4H  JJ*. I3//20X.3M  K*,I3) 

GO  TO  20 
30  WRITE  (6.7) 

WRITE (6*3) 

3 FORMAT (11  OH  PROGRAM  TO  COMPUTE  THE  FREE  SURFACE  DISTURBANCE  DUE  T1 
1 A STATIONARY  PULSATING  SOl'RCF  IN  WATER  OF  FINITE  DEPTH) 

WRITE (6*4) 

4 FORMAT (////✓) 

WRITE (6*5) 

5 FORMAT (22H  THE  INPUT  VALUES  ARE-) 

WRITE (6*6)  0*G*F*m*a*B*T 

6 FORMAT (23X*17H  SOURCE  STRENGTH*. F10.4//23X.22M  GRAVITATION  CONSTAN 
lTr,F10.4//23X,21H  PULSATING  FREQUENCY* *F10 .4//23X . 1 3H  WATER  DEPTH= 
2.P10.4//23X.22H  X LOCATION  OF  SOURCE=«F10.4//23X*22H  Y LOCATION  OF 
3 SOURCE*. FI 0.4//23X*6H  TIMP*,F10.4) 

WRITE (6* 1 1 ) STEPX.STEPY.STCPZ 

11  FORMAT (//23X.7H  STEPX* *F 1 0 .4//23X * 7H  STEPY**F10.4//23X.7H  STOPZ*.I 
13) 

WR I TE ( 6 » 7 ) 

7 FORMAT (1H1) 

WRITE(6*12)  SIGMA,GMU*ANGL1*SINF.C0EF1«C0EF1A,C0EF1B 

12  FORMAT  (7H  SIGMA* ,F  1 0 .4/5H  C-MU*.F10.6/7H  ANGL1  = *F10.2/6H  SINF**F10. 
16/7H  COEF1*,F10.4/8H  COEF1 A*,F10.4/BH  C0EF1B**F10.4) 

WRITE (6*8)  LAMBDA 

8 FORMATU8H  THF  WAVELENGTH  TS.1X.F10.4) 

WRITE (6*4) 

WRITE (6*9) 

9 FORMAT (5X *?H  X* 10X.2H  Y.11X.12H  WAVE  HEIGHT *4X*6H  LOOPS. 8X.5H  M (K ) 
1.10X.6H  COEF2.9X.5H  PHI1.9X.9H  PHI2A(K)) 

WRITE (6. 10) (X(II) .Y(II) .ETA(II) .NO(II) .PMU(II) .WAYOUT(II) .PHI1 (II) 
1 .PHI2A ( I I ) « I I»2*STOPZ ) 

10  FORMAT (lX.F10,4,2X,Fin«4,7X,F10»5.8X.I3.5X«F10»4,5X.F10»4.5X»F10»4 
1.6X.F10.6) 

20  STOP 
END 
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o r> 


SUBROUTINE  POOT1 (SIGMAH.FMt.N) 

N-0 

3 IF(SIGMAH-.7>  10*10*4 

4 IF (SIGMAH-1 .0)  13*60*5 

5 IF (SIGMAH-1 .2)  15*15*6 

6 FMU*SIGMAH 

7 GO  TO  21 

10  0*.5*SQRT(9.-12.*SIGMAH) 

Fmu*S0PT(1.5-D) 

GO  TO  ?1 
13  Fmu*1.19 
GO  TO  21 
15  FMU*1.2 
21  N=N*1 

A=FMU«(TANH(FMU) l-SIGMAH 

B=FM(J*  ( 1 . 0-TANH  (FMU)  *TANH (FMU) ) ♦TANH  (FMU) 

C*ABS(A/B) 

IF( (C/FMU)-. 00001)  50*50.30 
30  IF(N-IOO)  40.50*50 
40  FMU*FMU-A/B 
GO  TO  21 
50  RETURN 
60  FMU*1. 19967 
PFTUPN 
END 

SUBROUTINE  TO  COMPUTR  THE  REAL.  POSITIVE  ROOTS  OF  THE  EQUATION 

(M)*TAN(M)*P*0 

SUBOOUT I NE  POOT2 ( S I GM AH  * AMt . I ,K ) 

C TESTS  TO  SET  AN  INITIAL  VALUE  FOP  AMU 
IF  (K-2)  4.4,5 

4 AMU=( (SIGMAH*.01)/SIGMAH) *3. 1416/2. 

GO  TO  8 

5 G*K 

AMU=( (SIGMAH*.01)/SIGMAM) *3. 1416/2.4(6-2.0 >*3.1416 
GO  TO  R 

C START  OF  THE  NFWTON-RAPHSON  ITERATION 
R 1 = 0 
9 1*1*1 

T = AMU*T  AN ( AMU) *SIGMAH 

B= AMU* ( 1 .0*TAN ( AMU) *TAN ( AML ) ) *TAN ( AMU) 

R=ABS(T/B) 

C TEST  FOP  SI7E  OF  INCREMENTAL  ESTIMATE 
IF (P.LE.. 00001)  PFTUPN 
IF(I.GT.IOO)  RETURN 

amu=amu-t/b 

GO  TO  9 
ENO 
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Sample  Output 


MO«M«  TO  CO-.UTT  TMf  TUtT  SUOf.Cf  Ol'Tun.NCT  KM  To  . 1T.TIOM.Y  KJH.IINO  VXJ.CT  I.  ,.TT.  or  f | M I TT 


TMf  IN»"T  VALUCS  ••*- 

5OU0CC  ST»fH6TW«  1.0000 

c*Avrrorr<7««  c o*sto«t«  i*.*ooo 

POtSATIMG  r »f ouf Nf  »•  3.5000 

HM Tf»  OfPrum  A. oooo 
■ LOC0TIO*  Of  50U#Cf • 0.0000 

r Loano*  or  sou«cf»  -*.0000 


7 f**C  • 

*.O»90 

STCPI« 

1 .0000 

STCPr* 

0.0000 

STOP?* 

** 

5IGMA*  .IPO* 

r»*U»  ,30??57 

*«Gt I • I.)? 

5|NT«  ,TO?II0 

cocr 1 • .0*01 

cofrio-  i.oooo 

cocr 10*  .0*30 

TMt  «AVrv.FMGTM  IS  10. *0*0 


I 

T 

• *Vf  MfIGHT 

LOOPS 

M(R| 

cocr* 

PHJ» 

P*!*A«n 

0.0000 

0,0000 

.0 19?! 

5? 

*0.701* 

.000? 

.*1*1 

-.??1*0I 

l .0*00 

0.0000 

.01039 

1 3 

0 . ? 7 3 1 

.0*00 

.*101 

-.100?*? 

*.0"00 

0.0000 

.0150 7 

10 

* .0909 

.0150 

.*101 

-.1*70*0 

J.0000 

0.0000 

.01050 

0 

3.0*?9 

.0*59 

.*101 

-.004*55 

*.0000 

0.0000 

.0*301 

0 

3.0*79 

.0*59 

.*101 

-.05*199 

S. 0O00 

0.0000 

.013*0 

0 

3.0*79 

.0*59 

.*101 

-.0175  ?0 

*.0000 

0.0000 

• 00*7 1 

0 

3.0*79 

.0*59 

.*101 

-.0*5*1? 

?.o<*oo 

0.0000 

-.OOTTO 

0 

3.0*?9 

.0*59 

.*1*1 

-.01*9  1* 

0.0000 

0.0000 

-.01*70 

0 

3.0*79 

.0*59 

.*101 

-.011 1»5 

* . 00  00 

0.0000 

-.0??** 

0 

3.0*79 

.0*59 

.*101 

-.007051 

10. onoo 

0.0000 

- .0*5*0 

0 

3.0*79 

.0*59 

.*101 

- .0051*0 

1 1 .0000 

0.0000 

-.0**17 

0 

3.0*79 

.0*59 

.*1*1 

-.001*0* 

l?.0O00 

0.0000 

-.01919 

0 

3.0*79 

.0*59 

.*101 

-.00*1*9 

13.0000 

0.0000 

-.01131 

0 

3.0*79 

.0*59 

.*101 

-.001500 

1 * . 0O00 

0.0000 

-.00171 

7 

3.1*1* 

.051? 

.*1*1 

-.001055 

15.0000 

0.0000 

.0001? 

7 

3.1*1* 

.053? 

.*1*1 

-.000710 

i*.oooo 

0.0000 

.01005 

7 

3.1*1* 

.053? 

.*1*1 

-.000*70 

1T.0O00 

0.0000 

.0*30* 

7 

3.|?1* 

.0537 

.*1*1 

-.0001*1 

10.0000 

0.0000 

• 0/50? 

7 

3 . 1 ? 1 * 

.0537 

.*101 

-.000*10 

19.0000 

0.0000 

•0**77 

7 

i.m* 

.0517 

.*301 

-.0001*0 

?0.0000 

0.0000 

.0*00* 

7 

3.1*1* 

.0537 

.*1*1 

- .00009* 

*» .0000 

0.0000 

. 0 1 ?3* 

7 

3.1*1* 

.0517 

.*301 

-.000005 

**.0000 

0.0000 

• 00?* 1 

7 

3.1*1* 

.0537 

.*101 

- .0000** 

*1.©000 

0.0000 

-.00713 

7 

3.1*1* 

.0537 

.*1*1 

-.000010 

**.0000 

0.0000 

-.0100? 

7 

3.1*1* 

.0537 

.*  1*1 

-.0000*0 

*5.©*©© 

0.0000 

• • 0**53 

7 

3.1*1* 

.0537 

.*101 

-.00001* 

*0 . OflOO 

0.0000 

- • 0*509 

7 

3.1*1* 

.0517 

.*101 

-.000009 

/T.0O00 

0.00  JO 

-.0750* 

7 

3.1*1* 

.0537 

.*101 

-.000000 

*0.0000 

0.0000 

-.0*00? 

7 

3.1*1* 

.0517 

.*101 

-.00000* 

*9 . Oo  00 

0.0000 

-.01  1*1 

7 

3.1*1* 

.053? 

.*101 

-.000001 

30.0000 

0.0000 

-.00103 

7 

3.1*1* 

.053? 

• ? 10  1 

-.00000* 

31 .0000 

0.0000 

.0001* 

7 

3.1*1* 

.053? 

.*101 

-.000001 

3?.0000 

0.0000 

.01 5*0 

7 

3.1*1* 

.0537 

.*101 

-.000001 

31.0O00 

0.0000 

.0??01 

7 

3.1*1* 

.051? 

.*101 

-.000001 

U.OftOO 

0.0000 

• 0 *55* 

7 

3.1*1* 

.051? 

.*101 

-.000000 

35.0000 

0.0000 

. ®*5** 

7 

3.1*1* 

.051? 

.*101 

-.000000 

30.0000 

0.0000 

•0*1*7 

? 

3.1*1* 

.051? 

.*1*1 

-.000000 

3T.0000 

0.0000 

.01*10 

7 

3.1*1* 

.0517 

.*701 

-.000000 

30.0*00 

0.0000 

.00*0* 

7 

3.1*1* 

.051? 

.*101 

-.000000 

10. 0000 

0.0000 

-.00515 

7 

3.1*1* 

.051? 

.*101 

-.000000 

*0.0000 

0.0000 

-.01*)* 

f 

3.1*1* 

.053? 

.*103 

-.000000 
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